Process
- make buffer (400 m radius circles around well, which covers about 50
ha of land) for each well
source("./GitControlled/Codes/0_1_ls_packages.R")
source("./GitControlled/Codes/0_0_functions.R")
# /*===== point data for well location =====*/
data_w_LR_TB <- readRDS("./Shared/Data/WaterAnalysis/data_w_LR_TB.rds")
well_sf <- data_w_LR_TB %>%
unique(.,by="wellid") %>%
.[,.(wellid, longdd, latdd)] %>%
# --- the geographic coordinate in the well data is NAD83 (espg=4269)--- #
st_as_sf(., coords = c("longdd","latdd"), crs = 4269)
# /*===== make buffer around each well =====*/
well_buffer_sf <- well_sf %>%
#--- project to WGS UTM 14 (Cartesian 2D CS covering NE)---#
st_transform(32614) %>%
# --- make buffers--- #
# 400 m radius circles around well, which covers about 50 ha of land
st_buffer(., dist = 400)
ggplot()+
geom_sf(data=well_sf, size=0.5) +
geom_sf(data=well_buffer_sf, color="red", fill=NA)
well_buffer_sp <- as(well_buffer_sf, "Spatial")
- Then, apply
get_ssurgo_props() for each well
buffer.
get_ssurgo_props <- function(field, vars, summarize = FALSE) {
# field = as(filter(well_buffer_sf, wellid==112576), "Spatial")
# field = as(filter(well_buffer_sf, wellid==112591), "Spatial")
# field = as(filter(well_buffer_sf, wellid==231010), "Spatial")
# field= as(well_buffer_sf[5,], "Spatial")
# vars = soil_var
# summarize = TRUE
# Get SSURGO mukeys for polygon intersection
ssurgo_geom <-
SDA_spatialQuery(
geom = field,
what = "geom",
db = "SSURGO",
geomIntersection = TRUE
) %>%
st_as_sf() %>%
dplyr::mutate(
# area = as.numeric(st_area(.)), #this causes an error for some polygons
area = .$area_ac,
area_weight = area / sum(area)
)
# --- check --- #
# ggplot()+geom_sf(data=ssurgo_geom, aes(fill=factor(mukey)))
# Looks like, in one buffer (surrounding a well), there are multiple polygon data where each of them are associated with a unique mukey
# `area_weight` indicates a portion of each polygon to entire butter
# Get soil properties for each mukey
mukeydata <-
get_SDA_property(
property = vars,
method = "Weighted Average",
mukeys = ssurgo_geom$mukey,
top_depth = 0,
bottom_depth = 150
)
ssurgo_data <- left_join(ssurgo_geom, mukeydata, by = "mukey")
if (summarize == TRUE) {
ssurgo_data_sum <-
ssurgo_data %>%
data.table() %>%
.[,
lapply(.SD, weighted.mean, w = area_weight, na.rm = TRUE),
.SDcols = vars
]
return(ssurgo_data_sum)
} else {
return(ssurgo_data)
}
}
Example: SSURGO
polygon data for an example field
library(soilDB)
library(aqp)
library(sp)
ex_well_bf <- well_buffer_sf[5,]
ex_well_bf_sp <- as(ex_well_bf, "Spatial")
ssurgo_geom <-
SDA_spatialQuery(
geom = ex_well_bf_sp,
what = "geom",
db = "SSURGO",
geomIntersection = TRUE
) %>%
st_as_sf()%>%
dplyr::mutate(
area = as.numeric(st_area(.)),
area_weight = area / sum(area)
)
# /*===== visualization =====*/
ggplot()+
geom_sf(data=ssurgo_geom, aes(fill=factor(mukey)), color="blue", alpha=0.5)+
geom_sf(data= well_sf[5,], fill=NA)
Download SSURGO
data
- This is the code I used to download SSURGO data for my data.
# === Target variables === #
soil_var <- c(
# --- "Horizon" --- #
"sandtotal_r", # a fraction of sand (0.05mm to 2.0mm)(%)
"claytotal_r", # a fraction of clay (%)
"silttotal_r", # a faction of silt (%)
"ksat_r", # hydraulic conductivity (um/m)
"awc_r", # available water capacity (cm/cm)
# --- "Component" --- #
"slope_r" # The difference in elevation between two points, expressed as a percentage of the distance between those points. (SSM), "component"
)
# each filed (buffer) data is passed on get_ssurgo_props()
res_ssurgo <- lapply(
seq_len(nrow(well_buffer_sp)),
function(x) {
print(paste0("working on wellid: ", well_buffer_sp[x, ]$wellid, " field"))
# x=1
get_ssurgo_props(
field = well_buffer_sp[x,],
vars = soil_var,
summarize = TRUE
) %>%
.[, wellid := well_buffer_sp[x, ]$wellid]
}) %>%
bind_rows()